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Introduction 

The  objective  of  this  project  is  to  develop  an  automatic  oceano¬ 
graphic  image  interpretation  system.  This  system  will  aid  the  operator 
in  identifying  various  oceanic  features  and  will  predict  the  dynamic 
behavior  of  these  features.  Presently  these  tasks  are  handled  by  experts 
in  oceanography  and  obviously,  a  time  consuming  process.  The  pro¬ 
posed  project  has  been  divided  into  three  phases  :  In  Phase  I  the 
features  present  in  the  oceanographic  image  are  identified.  In  this  phase 
an  efficient  feature  labeling  algorithm  will  be  implemented.  In  Phase  II 
an  Expert  System  that  helps  in  predicting  the  dynamic  behavior  of  the 
oceanic  features  will  be  developed.  In  Phase  III  these  two  modules  will 
be  integrated  and  the  system  will  be  tested  with  real-time  data. 

Current  Progress 

The  design  and  development  of  Phase  I  started  in  May  1988.  Vari¬ 
ous  techniques  for  feature  labeling  were  studied.  As  a  result  of  a  series 
of  discussion  with  Ron  Holyer  and  Matthew  Lybanon  ,  it  was  decided 
to  use  non-linear  probabilistic  relaxation  labeling  scheme  for  the 
oceanographic  feature  labeling  problem.  We  believe  that  this  approach 

V 

is  the  first  of  its  kind  ,  namely  ,  the  application  of  relaxation  labeling 
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for  identifying  oceanic  features.  The  design  of  the  feature  labeling 
module  was  completed  in  June  1988.  The  implementation  of  the  feature 
labeling  algorithm  was  initiated  in  June  1988.  As  a  first  step,  the  a 
priori  probabilities  of  all  the  features  are  estimated  and  tested  on  a 
reference  image.  The  results  obtained  in  the  first  step  are  published  in 
5th  International  conference  on  Intelligent  Information  Processing  Sys¬ 
tems,  to  be  held  in  Jan  1989.  The  second  step  in  the  relaxation  scheme 
called,  "iterative  updating"  is  ready  for  implementation.  The  design  of 
this  sub-module  was  completed  in  August  1988.  It  is  proposed  to  com¬ 
plete  this  implementation  by  October  1988.  The  experimental  results 
were  published  in  two  parts  in  two  conferences.  For  more  information 
on  these  results  please  refer  : 

1.  A  Technique  for  Feature  Labeling  in  Infrared  Oceanographic 
Images  ,  appeared  in  Fifth  International  Conference  on  IIPS,  Jan 
1989. 

2.  An  Expert  System  for  Interpreting  Mesoscale  Features  in 
Oceanographic  Satellite  Images,  appeared  in  SPIE  Conference  on 
Applied  Artificial  Intelligence,  March  1989. 

The  design  and  development  of  the  Phase  II  of  the  project  was  ini¬ 
tiated  in  December  1988.  The  ground  work  for  the  design  of  Phase  II 
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has  already  been  started.  We  have  started  analyzing  the  various  com¬ 
ponents  involved  in  an  Expert  System,  especially  for  the  oceanographic 
application.  The  objective  of  the  Phase  II  is  to  develop  a  powerful 
expert  system  that  incorporates  the  nature  of  the  oceanic  features  as 
facts  and  rules.  The  expert  system  will  make  use  of  the  results  obtained 
from  Phase  I.  The  two  modules  will  then  be  integrated  in  Phase  III. 


A  Report  on  Feature  Labeling  Software 

The  software  modules  are  written  in  C.  In  this  report,  we  give  some  general 
comments  about  the  program  and  various  assumptions  made.  The  software 
modules  are  to  run  on  VAX-VMS  system.  Several  system-related  functions 
are  used  and  hence  modifications  are  required  if  necessary  to  port  to  another 
machine.  The  modules  are  written  in  standard  C.  The  following  paragraph 
describes  the  source  files  and  input  files  used  in  the  software: 

1.  Al.C 

Source  file  containing  the  first  part  of  the  a  priori  program.  It  cal¬ 
culates  the  mean,  covariance  vectors  etc 
INPUT  :  there  are  FIVE  input  files  for  this  module. 

a)  IMAGE.DAT  -  the  actual  image  (512  X  512)  in  the  line 
oriented  format.  To  prepare  this  image  take  the  (512  X  512)  I  S 
file;  give  it  as  input  to  the  C’DISKTRANSFER  function;  rename 
the  output  as  IMAGE.DAT. 

b)  CLUSTER.DAT-the  output  of  the  cluster  shade  edge  detector 
algorithm.  Again  this  file  HAS  TO  BE  i  the  line  oriented  format. 

c)  ANGLE.DAT  -  this  file  contains  the  angle  subtended  by  the 
vector  (X,Y)  to  the  X-axis.  This  file  is  prepared  only  once  and 
hence  need  to  be  just  copied  while  executing  A 1  .EXE  file. 

d)  DIST.DAT  -  this  file  contains  the  magnitude  of  the  vector 
(X,Y).  this  is  very  similar  to  ANGLE.DAT  ;  need  to  be  just 
copied  while  executing  Al.EXE  file. 
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e)  OBJ.DAT  -  this  file  contains  the  information  about  the 
objects  that  are  to  be  labeled  in  the  image.  This  is  a  structured 
file.  The  information  about  the  objects  are  got  from  an  external 
data  file.  To  prepare  this  file,  a  program  called  RDATA.EXE 
has  to  be  executed  for  every  new  image  and  data  set.  The  details 
about  this  is  given  laier  in  thius  report 

OUTPUT  :  a  file  called  MEAN.DAT  which  contains  the  mean 
and  covariance  vectors  for  the  given  image  and  data. 


2.  A2.C 


Source  file  containing  the  second  part  of  the  a  priori  program. 

INPUT  :  There  are  SIX  input  files.  They  are  IMAGE.DAT, 
CLUSTER.DAT,  ANGLE.DAT,  DIST.DAT,  OBJ.DAT  and 
MEAN.DAT. 

OUTPUT:  The  output  of  this  module  is  the  a  priori  image.  It  is 
stored  in  a  file  called  APR.IMAGE.  This  file  will  be  used  subse¬ 
quently  in  the  labeling  process. 

3.  LABEL.C 

Source  file  containing  the  last  part  of  the  labeling  program, 
namely,  assigning  labels  to  the  objects. 

INPUT:  -  TWO  input  files  :  CLUSTER.DAT  and  APR.IMAGE 
(prepared  by  A2.EXE). 
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OUTPUT:  -  A  file  called  LABEL.IMAGE  which  is  the  final  out¬ 
put  of  the  feature  labeling  software  and  it  is  NOT  DISPLAY- 

ABLE.  The  output  image  is  still  in  the  line  oriented  format.  The 

2 

last  step  is  to  convert  this  image  into  I  S  file  format  ,  use 
C’DISKENTER  function.  We  can  call  the  displayable  output  as 
LABEL.DIS. 

4.  RDATA.C 

Source  file  that  converts  the  external  data  file  (containing  the  points  on  the 
GULF  STREAM  etc.)  to  the  format  understandable  by  the  feature  labeling 
module. 

INPUT  :  -  an  EXTERNAL  DATA  file  (eg.,  mpr8605) 

OUTPUT:  -  OBJ.DAT  which  is  used  by  A1  and  A2  programs.  If 
you  feel  that  there  are  more  objects  than  needed,  you  can  edit  the 
OBJ.DAT  file  (which  is  in  the  text  format)  with  the  help  of  the 
information  given  below  as  a  note. 

NOTE  ON  THE  STRUCTURE  OF  OBJ.DAT: 

The  objects  to  be  labeled  are  of  two  types:  one  with  all  the  pixels  on  it  listed 
(eg.,  GULF  STREAM,  SHELL  etc)  and  the  other  with  center  and  radius 
given  (eg.,  WCR,  CCR).  The  first  type  is  given  a  code  0  and  the  second  type  a 
code  1  in  the  OBJ.DAT  file  structure.  Hence  the  structure  is  :  code  of  the 
object  followed  by  the  data.  In  the  case  of  object  type  0,  the  data  contains  all 
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the  pixel  coordinates  and  in  the  case  of  object  type  1,  the  data  contains  the 
center  and  radius. 


HOW  TO  RUN  THE  PROGRAM  ? 

First  compile  these  source  files  using  CC  command  and  use  LINK  command 
to  prepare  the  EXE  files.  For  example  to  create  Al.EXE  use  the  following 
sequence  of  commands. 

CC  Al.C 
LINK  A 1 

Run  the  programs  Al.EXE,  A2.EXE  and  LABEL.EXE  in  the  correct 
sequence  to  get  the  final  output. 


For  further  information  please  contact  Prof.S.S.Iyengar  or  N.Krishnakumar  at 
LSU  (Ph#  504-388-1249). 
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ABSTRACT 

Thermal  infrared  images  of  the  ocean  obtained  from  satellite  sensors  are  widely  used  for  ±e  study  of  ocean 
dynamics.  The  derivation  of  mesoscale  ocean  information  from  satellite  data  depends  to  a  large  extent  on  the 
correct  interpretation  of  infrared  oceanographic  images.  The  difficulty  of  the  image  analysis  and  understanding 
problem  for  oceanographic  images  is  due  in  large  part  to  the  lack  of  precise  mathematical  descriptions  of  the 
ocean  features,  coupled  with  the  time  varying  nature  of  these  features  and  the  complication  that  the  view  of  the 
ocean  surface  is  typically  obscured  by  clouds,  sometimes  almost  completely.  Towards  this  objective,  the  present 
paper  describes  a  hybrid  technique  that  utilizes  a  nonlinear  probabilistic  relaxation  method  and  an  expen  system 
for  the  oceanographic  image  interpretation  problem.  A  unified  mathematical  framework  that  helps  in  solving  the 
problem  is  presented.  This  paper  highlights  the  advantages  of  using  the  contextual  information  in  the  feature 
labeling  algorithm.  The  paper  emphasizes  the  need  for  the  feedback  from  the  high  level  modules  to  the  intermedi¬ 
ate  modules  in  an  automatic  image  interpretation  system.  The  paper  presents  some  important  results  of  the  series 
of  experiments  conducted  at  Remote  Sensing  Branch,  NORDA,  on  the  NOAA  AVHRR  imagery  data. 

Key  words:  feature  labeling,  feature  extraction,  oceanic  features,  edge  detection,  knowledge  based  systems, 
expert  system,  relaxation,  infrared  imagery. 


1,  INTRODUCTION 

Satellite-bome  sensors  potentially  offer  many  advantages  for  the  study  of  oceanic  processes.  They  provide 
global  synoptic  measurements  of  various  oceanic  surface  properties,  in  contrast  to  the  local  measurements,  possi¬ 
bly  at  a  range  of  depths,  provided  by  conventional  oceanographic  measurement  techniques.  Thermal  infrared  (IR) 
images  of  the  ocean  obtained  from  satellite  sensors  are  widely  used  for  the  study  of  ocean  dynamics.  Fig.  1  shows 
a  sample  infrared  image  of  the  Gulf-Stream  obtained  from  the  Advanced  Very  High  Resolution  Radiometer 
(AVHRR)  aboard  the  NOAA-7  satellite.  Brightness  in  this  infrared  image  is  inversely  proportional  to  the  ocean 
surface  temperature  (dark  areas  represent  warmer  temperatures  and  light  areas  represent  colder  temperatures). 
Vortices  (areas  of  closed  circulation)  within  this  turbulent  flow  pattern  are  called  eddies.  The  Gulf  Stream  and  its 
associated  eddies  are  examples  of  mesoscale  features  ("mesoscale"  is  the  name  commonly  applied  to  the  features 
existing  on  spatial  scales  of  the  order  of  50  to  300  km).  Mesoscale  features  are  important  to  the  study  of  ocean 
dynamics,  to  the  fisheries  and  to  many  other  diverse  interests. 
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C.irer.’  -  ue  ar.a’.y ms  techniques  rely  on  human  interpretation  of  the  satellite  imagery.  Human  interpreta- 
>  ..  varied  m  its  level  of  expertise  ar.d  is  highly  labor-intensive.  With  the  proliferation  of  high 

Ad-. .ittcea  Very  High  Resolution  Radiometer  image  applications  .  u  becomes  highiy  desirafcie  for  certain 
.  to  move  tram,  the  labor-intensive  manual  interpretation  of  infrared  imagery  towards  a  capability  for 
mated  interpretation  of  these  images  The  complete  automation  of  the  oceanographic  image  interpretation 
is  rmr.oo'o.  not  fea.able.  but  one  can  begin  to  address  certain  subsets  of  the  probiem  with  the  present-, Tv 
ge  processing  and  artmciai  intelligence  techniques.  This  was  the  motivation  for  the  work  reported  by  Lyb-ner. 
h.  •  1 1.  in  the  development  of  a  prototy  pe  oceanographic  expert  system. 
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.".era!  prev  :ous  studies  have  addressed  the  automation  of  the  analysis  of  infrared  satellite  imagery  for 
...'e  features  Gerscn  and  Gabrc.sk:  [2]  and  Gerson  et  ai . [ 3 1  investigated  the  detection  of  the  Gulf  Stream,  it 
images  from  the  Geostationary  Operational  Environmental  Satellite  (GOES).  Gerson  ar.d  Gabrcsk;  [dl 
uerarchicaJ  approach  where  16.xi6  pi.xei  (IZSxiZS  km)  "frames"  within  the  image  were  evaluated  for  the 
:tv  o:  c,  r.taimng  the  Gulf  Stream.  Frames  flagged  as  Gulf  Stream  possibilities  were  then  further  evaluated 
mine  the  exact  location  of  the  Stream  within  the  frame  by  looking  at  statistics  based  on  5.0  pixel  "local 
: meeds'.  As  an  outgrowth  of  the  work  reported  in  (2!  and  [3),  Coulter  [-]  performed  automated  feature 
1  studies  using  the  higher  resolution  Advanced  Very  High  Resolution  Radiometer  data,  ianowicz  [5]  stu- 
lutcmatic  detection  of  the  Gulf  Stream  eddies  using  Advanced  Very  High  Resolution  Radiometer  data, 
i  used  a  region  adjacency  graph  to  define  spanai  relanonslups  between  elementary  connected  regions  of 
crav  level  called  atoms.  Althoucn  sansfactorv  emulation  of  human  extraction  of  eddv  structure  is 


this  method,  N ichci (6 j  did  point  cut  that  net  ail  enclosed  uniform  areas  identified 


A  t. 
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end  to  real  ocean  structure.  The  authors  [7]  nave  approached  this  preoiem  and  shown  some  preliminary 
cn  the  relaxation  based  feature  labeling  scheme,  the  development  of  a  prototype  expert  system  and  axe 
:e  between  these  two  modules. 


2.  MOTIVATION'  OF  THE  PRESENT  WORK 


In  this  section  we  present  the  motivation  behind  the  present  work.  Our  primary  objective  is  to  build  3 
powerful  automatic  image  interpretation  system  for  oceanographic  satellite  images.  In  order  to  make  this  difficult 
probiem  tractable,  we  divide  the  problem  into  two  parts:  feature  labeling  problem  and  the  development  of  an 
expert  system.  It  is  ciear  that  the  performance  of  the  labeling  algorithm  depends  heavily  on  the  low  level  image 
processing  algorithms.  Particularly,  the  output  of  an  edge  detector  algorithm  plays  a  major  role  in  the  feature 
labeling  process  In  view  of  this,  a  new  efficient  edge  detector  algorithm  proposed  by  Holyer  and  Peckinpaugh  [8] 
is  employed  in  our  feature  labeling  technique.  It  is  known  that  the  conventional  edge  derivative  operators  are  very 
sensinve  to  noise  and  are  not  suitable  for  analyzing  oceanographic  satellite  image  The  new  edge  detector  algo¬ 
rithm  proposed  by  Holyer  and  Peckinpaugh  [8]  is  based  on  the  gray  level  co-occurrence  matrix,  which  is  com¬ 
monly  used  in  image  texture  analysis.  This  algorithm  [8]  is  found  to  exhibit  the  characteristics  of  fine  structure 
rejection  while  retaining  edge  sharpness.  In  this  paper  we  focus  on  the  preliminary  results  of  a  feature  labeling 
algorithm  and  the  development  of  an  expert  system  to  aid  the  interpretation  of  these  features.  In  the  follownng 
paragraph  we  briefly  outline  the  various  steps  taken  to  solve  the  second  subproblem. 

The  objective  of  an  oceanographic  expert  system  is  to  correctly  interpret  the  dynamics  of  the  ocean  process 
with  minimal  human  interacnon.  Towards  this  objective,  the  development  of  a  powerful  oceanographic  expert 
system  is  under  way  at  Remote  Sensing  Branch,  NORDA.  We  have  already  developed  a  prototype  expert  system 
for  Gulf  Stream  regional  dynamics  [1].  The  current  activities  focus  on  building  an  expert  system  that  makes  use 
of  the  information  about  the  position  of  mesoscale  features  obtained  from  the  relaxation  labeling  scheme.  One 
element  of  an  expert  system  is  a  database  of  knowledge  about  the  subject  matter,  with  the  knowledge  represented 
in  a  form  suitable  for  manipulation  by  the  "inference  engine"  of  the  expert  system  (i.e.,  the  logical  and  heuristic 
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procedures  for  solving  problems  in  the  problem  domain).  With  the  help  of  the  knowledge  gained  from  discussions 
at  NORDA  and  from  the  literature  in  oceanography,  a  knowledge  base  about  the  nature  of  the  mesoscale  features 
(occurrence,  mean  lifetime,  movement,  etc.)  was  built  [1].  Presently,  we  are  in  the  process  of  designing  an  expert 
system  which  makes  use  of  positional  information  about  the  features(obtained  from  the  relaxation  labeling)  and 
the  know  ledge  base  of  mesoscale  information. 

This  paper  presents  a  relaxation  labeling  scheme  to  label  features  in  thermal  infrared  images  obtained  from 
satellite.  The  technique  presented  m  this  paper  exploits  the  advantages  of  using  contextual  information  in  the 
labeling  algorithm.  The  remainder  of  the  paper  is  organized  as  follow's:  Section  3.0  gives  an  overview  of  the 
architecture  of  the  proposed  image  interpretation  system.  Section  4.0  introduces  the  probabilistic  relaxation 
scheme  and  briefly  discusses  the  application  of  this  scheme  to  the  oceanographic  labeling  problem.  In  section  5.0 
we  develop  a  mathematical  framework  which  is  necessary  for  our  labeling  problem.  Section  6.0  discusses  the  per¬ 
formance  of  cur  technique  and  describes  the  steps  that  are  involved  in  implementing  this  technique.  Section  7.0 
describes  the  prototype  expert  system  and  the  interface  between  the  expert  system  and  the  feature  labeling 
module,  section  S.O  concludes  with  all  the  important  features  of  this  technique  and  future  extensions  to  the  avail¬ 
able  scheme. 


3.  ARCHITECTURE  OF  THE  INTERPRETER 

In  this  section,  we  briefly  describe  the  various  components  involved  in  the  image  interpretation  system.  The 
block  diagram  shown  in  figure. 2  depicts  the  underlying  architecture  cf  the  proposed  system.  The  IR  image 
obtained  from  the  satellite  is  the  input  to  the  feature  labeling  module.  The  feature  labeling  module  consists  of  two 
submodules  :  one  estimating  the  initial  feature  probabilities  and  one  implementing  the  iterative  updating  scheme. 
The  initial  pixel  probabilities  are  assigned  with  the  help  of  the  previous  analysis  and  the  ground  truth  data.  These 
probability  values  are  then  iteratively  updated  using  the  nonlinear  probabilistic  relaxation  scheme.  Finally,  the 
relevant  features  that  are  present  in  the  image  are  labeled  with  the  aid  of  a  decision  process. 

The  labeled  features  are  fed  to  the  expert  system  module.  This  module  is  divided  into  two  submodules.  The 
first  submodule  consists  of  a  knowledge  base  and  an  i  Terence  engine.  Tne  knowledge  base  has  a  thorough  collec¬ 
tion  of  facts  and  rules  about  the  mesoscale  features  svch  as  gulf  stream,  warm  core  rings  etc.,  Based  on  the  posi¬ 
tional  information  obtained  from  the  feature  labeling  module  and  the  knowledge  base,  the  inference  engine 
attempts  to  interpret  the  dynamics  of  these  features.  E  asically  the  inference  engine  is  a  pattern  matching  module. 
The  second  submodule  evaluates  the  consistency  of  t  ;e  labeling  process  with  the  help  of  the  knowledge  base  as 
well  as  from  the  input  from  the  operator.  This  evaluation  results  in  assigning  a  confidence  factor  to  each  of  the 
features  detected.  This  factor  is  then  "fed  back"  to  the  feature  labeling  module,  in  addition  to  the  previous  data 
analysis,  to  improve  the  consistency  of  the  labeling  process  in  the  Mature. 


4.  RELAXATION  PROCESS 

An  important  research  area  in  image  analysis  and  image  interpretation  technology  is  the  development  of 
methods  that  blend  contextual  information  with  conventional  image  processing  algorithms.  A  literature  survey 
clearly  indicates  that  such  a  hybrid  approach  yields  good  results.  Relaxation  labeling  is  one  such  process.  Relax¬ 
ation  labeling  has  been  applied  to  a  variety  of  image  processing  problems  e.g.,  linear  feature  enhancement  [9], 
edge  enhancement  [10],  image  enhancement  [11],  pixel  classification  [12,13].  A  recent  survey  article  by  Kitten 
and  Illingworth  [14]  on  relaxation  labeling  highlights  the  importance  of  this  area  of  research.  The  survey  [14]  also 
points  out  the  advantages  and  possible  applications  of  relaxation  methods.  More  importantly,  the  relaxation  label¬ 
ing  approach  was  elegantly  described  by  Rosenfeld  et  al.  [15]  who  investigated  the  problem  of  labeling  the  sides 
of  a  triangle  and  proposed  a  set  of  schemes  to  solve  the  problem.  The  paper  [15]  concluded  with  the  result  that  the 
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r.or.lmear  probabilistic  relaxation  schemes  yield  better  results  than  the  others.  The  labeling  algorithm  presented  in 
this  paper  is  based  on  the  non  linear  probabilistic  relaxation  technique. 

The  goal  of  the  relaxation  process  is  to  reduce  the  uncertainty  (and  improve  the  consistency)  in  the  assign¬ 
ment  of  one  of  the  labels  to  each  object  in  a  set  of  related  objects.  Ln  the  oceanographic  feature  classification  prob¬ 
lem,  the  classes  are  the  various  oceanographic  features,  namely  the  north  and  south  walls  of  the  Gulf  Stream,  cold 
eddy,  wam:  eddy,  shelf  front  and  coastal  boundary.  Refer  to  figure. 3  to  identify  the  positions  of  these  oceanic 
features  in  a  typical  image.  The  objects  are  the  individual  pixels  in  a  set  of  registered  multi-temporal  images. 
The  uncertainty  could  be  due  to  the  cloud  cover  or  the  overlap  of  the  features,  features  not  belonging  to  one  of  the 
classes,  noise  Ln  the  image,  or  other  factors.  In  this  paper,  we  are  attempting  to  label  mesoscale  features,  but  the 
ccean  exhibits  variability  on  all  spatial  scales.  Thermal  structure  on  scales  smaller  than  mesoscale  will  interfere 
with  the  mesoscale  fearure  labeling  process.  The  underlying  mathematical  framework  necessary  for  the  relaxation 
labeling  method  is  described  in  the  next  section  of  the  paper. 


5,  MATHEMATICAL  FRAMEWORK 

Let  A  =  }  be  the  set  of  possible  labels  that  may  be  assigned  to  each  pixel  x  in  the  IR  image.  Also  we 

let  px(x,r)  denote  the  probability  that  the  pixel  at  x(ij)  and  time  instant  t,  belongs  to  the  object  X  after  k.  itera¬ 
tions  of  the  relaxation  algorithm.  Note  that  the  probabilities  are  functions  of  time  unlike  the  conventional  pixel 
relaxation  labeling  schemes  where  the  probability  is  a  function  of  position  alone.  This  allows  the  relaxation  label¬ 
ing  algorithm  to  utilize  temporal  continuity  to  reduce  the  ambiguity  in  labeling.  The  ambiguity  may  arise  due  to 
noise  (cloud  cover,  for  example). 

There  are  two  steps  in  executing  the  probabilistic  relaxation  algorithm.  In  the  first  step,  a  priori  probabilities 
are  evaluated  with  the  help  of  ground  truth  data  and  /  or  a  previous  but  recent  mesoscale  analysis.  In  the  second 
step,  these  a  priori  probabilities  are  iteratively  updated  (relaxation)  until  a  consistent  labeling  is  reached.  We  now 
discuss  these  two  steps  in  detail. 

Step  1:  Estimating  the  a  priori  probabilities 

Let  px  ,0  denote  the  a  priori  value,  that  is,  the  probability  that  pixel  x(ij)  at  time  t  belongs  to  the  object  X  at 
the  zeroth  iteration.  The  Bayesian  probability  equation  is  used  to  evaluate  this  value.  The  equation  (5.1)  is  used  to 
calculate  p\{x,t). 


P\(x,  0 


p(x,t\X)  P(k) 
Zx  p(x,t  |X)  P(k) 


(5.1) 


where  p  (x,t  i  denotes  the  conditional  density  function  and  P(X)  the  probability  of  occurrence  of  the  object  X. 

To  evaluate  'he  ccc-Ation al  density  function  p(x,t  |  X),  a  set  of  parameters  is  measured  at  the  pixel  x(ij).  Let  X 
denote  the  par  uneter  vector.  The  following  parameters  are  used  to  form  the  vector  X: 

(1)  vector  from  origrn  to  pixel  x(ij),  both  the  magnitude  and  direction. 

(2)  gray  scale  intensity  value  at  the  pixel  x(ij). 
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(3)  the  edge  magnitude  (Section  6. 1  presents  the  chosen  edge  operator  algorithm). 


For  each  object,  the  mean  vector  fix  and  the  covariance  matrix  Zx  are  computed.  Also  it  is  assumed  that  the  condi¬ 
tional  densicy  function  follows  a  normal  distribution.  Hence  the  conditional  density  function  p(x,t\h )  is 
evaluated  using  equation  (5.2). 


P(x,t  |  X)  =  (2ji |  Zx | )  *ex; 


2 


(x-pxy^c'c*-^) 


(5-2) 


To  compute  P  (X)  ,  relative  areas  of  the  objects  are  considered.  The  number  of  pixels  in  the  object  X  is  «x-  Then 
P  (X)  can  be  calculated  using  equation  (5.3). 


P(X)  = 


£x*x 


(5.3) 


Step  2:  Iterative  updating  algorithm 

We  now  discuss  the  probability  updating  rule.  The  new  estimate  of  the  probability  of  X  at  x(i  j)  is  given  by  (5.4). 


pi+l(x,t)  =  ax 


pt(x,r)(l-H7£(*)) 

(!+?£(*)) 


+  (1-a x)/>x(*/) 


(5.4) 


where  q  J  (x )  is  called  the  update  factor  and  ax  is  called  the  temporal  weighting  function. 

The  method  of  estimating  these  factors  is  illustrated  below.  The  updating  factor  for  the  estimate  px(x,r)  at  the 
fcth  iteration  is  given  by  equation  (5.5). 


qi{x)  =  ^-Iylxr>xix,y)pi{y)  (5-5) 

til 


where  m  is  the  number  of  objects.  In  this  equation,  denote  compatibility  coefficients.  These  coefficients 

are  computed  as  in  [13.15J.  According  to  the  relaxation  scheme,  r^xy)  is  a  measure  of  the  probabilistic  compa¬ 
tibility  between  label  X  on  point  x  and  label  X'  on  pointy,  and  has  the  following  characteristics: 
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I 

(1)  If  X  on  x  frequently  co-occurs  with  X  on  y,  then  rxx.ix  ,y )  >  0  ,  and  if  they  always  co-occur,  then  rxx<*  ,y )  « 
1. 

(2)  If  A.  on  x  rarely  co-occurs  with  X  ony,  then  rxx<x,y )  <  0  ,  and  if  they  never  co-occur,  then  r^{x,y)  -  -  I. 

(3)  If  X  on  x occurs  independently  of  X  ony,  then  rxx<x,y)  =  0, 

These  compatibility  coefficients  are  computed  using  equation  (5.6). 


fxxU.y) 


cov  (X,X  ) 
o(X)  o(X') 


(5.6) 


where  cov  (A^X )  denotes  the  covariance  of  two  events,  namely,  the  pixel  x  is  assigned  the  label  X  and  the  pixel  y  is 
assigned  the  label  X ,  ahd  is  given  as: 


cov  (X,X )  =  p  (k&  X  )  -  p  (X)p  (X ) 


(5.7) 


The  joint  probability  p  (KA  X )  can  be  calculated  using  the  position  vectors  of  the  pixels  x  and  y.  Let  m  and  m’  be 
the  mean  position  vectors  of  the  objects  X  and  X  respectively.  Then  the  joint  probability  can  be  computed  using 
simple  functions  as  given  by  equations  (5.8a)  and  (5.8b). 


p  (X  &  X')  = 

if  X  *  X 

(5.8a) 

m  —  m 

P(k  a')  =  i  - 

x+m 

if  X  =  X" 

(5.8b) 

« 

The  new  estimate  of  the  probability  of  X  at  x(i  j)  is  : 


P  £+l  C*  .0  =  «x 


ffcfr.O  (!+?£(*)) 
£xp£(*,0  (!+?£(*)) 


+  (l-a\)Px(xS) 


(5.9) 


where  p\(x  / )  denotes  the  probability  value  at  the  pixel  x(ij)  at  time  instant  t’,t’  <  t  and  ax  denotes  the  temporal 
weighting  function  for  the  object  X.  The  function  ax  determines  the  weight  that  is  given  to  the  current  probability 
values  and  the  factor  (1— ax)  associates  a  weight  to  the  probability  values  calculated  at  previous  time  instant  t’ . 

To  estimate  the  new  probability  value,  either  (a)  the  probability  value  at  the  previous  time  frame  t’  or  (b)  a 
set  of  probability  values  at  time  frames  t ,  t” ,  t'"  etc.,  can  be  used.  In  the  present  analysis  and  the 
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implementation,  method  (a)  is  chosen.  Also  it  is  proposed  to  implement  method  (b)  and  compare  improvements 
(if  any)  in  the  feature  labeling.  The  temporal  weighting  function  can  be  determined  with  the  help  o*the  previous 
analysis  of  the  images.  It  is  understood  that  this  temporal  weighting  function  can  not  be  constant  in  the  real  time 
situation.  For  example  in  some  cases,  the  performance  of  the  labeling  may  be  less  satisfactory  and  correspond¬ 
ingly  the  temporal  weighting  function  should  be  given  a  relatively  smaller  value.  Hence  the  selection  of  this  func- 
uon  itself  may  be  a  problem  in  some  cases.  One  may  use  different  values  for  this  function  and  the  performance  of 
the  labeling  can  be  compared.  The  iterative  updating  is  terminated  when  the  difference  between  the  probability 
values  at  the  k  th  and  the  (k  +  1)  th  iterations  is  very  small  (say  less  than  0.1  %). 


6.  DISCUSSION  OF  THE  TECHNIQUE 

The  implementation  of  the  above  mentioned  technique  is  carried  out  in  two  stages. 

Stage  I 

At  the  first  stage,  the  a  priori  probabilities  are  estimated  using  a  manually  prepared  mesoscale  analysis  from 
a  time  period  of  five  days  prior  to  the  test  image. 

Tne  objects  present  in  the  oceanographic  IR  image  are  identified  with  the  help  of  the  previous  analysis.  The 
shapes  (or  boundaries)  of  these  objects  are  then  determined  and  represented  by  a  regular  polygon.  The  parameter 
values  are  then  computed  for  all  the  pixels  inside  these  polygons.  A  reasonably  accurate  algorithm  which  finds 
out  whether  a  pixel  is  inside  a  polygon  or  not  is  implemented. 

From  the  parameter  vector  X  for  a  pixel  x(i  j),  the  mean  vector  p*  is  calculated.  Also  for  each  object  X  ,  the 
covariance  matrix  Z^  is  computed.  Equation  (5.2)  is  used  to  compute  the  conditional  density  function p{x,t  |X). 
Finally  the  initial  probability  p£(x,t)  is  computed  using  the  equation  (5.1).  The  output  of  the  first  stage  is  shown 
in  figure.4.  Fig.4  shows  the  output  of  the  labeling  module  with  all  five  features  labeled.  Refer  to  figure.3  to  com¬ 
pare  the  labeling  with  the  hand  segmented  image.  To  illustrate  the  labeling  process  five  different  features  are  con¬ 
sidered.  They  are:  two  parts  of  the  north  wall,  south  wall  and  two  warm  eddies.  The  labeling  of  these  five  features 
is  shown  in  figure  4.  The  individual  pixels  in  the  oceanographic  image  are  labeled  depending  on  the  initial  proba¬ 
bility  values  obtained  in  the  first  stage. 

6.1  DISCUSSION  OF  CLUSTER  SHADE  EDGE  ALGORITHM 

In  this  section,  we  discuss  the  features  of  the  edge  detection  algorithm  proposed  by  Holyer  and  Peckinpaugh 
[8],  The  motivation  behind  the  development  of  such  a  new  edge  detector  algorithm  is  to  aid  the  analysis  of 
oceanographic  satellite  images.  The  popular  derivative-based  edge  operators  viz.  Sobers  operator  are  shown  to  be 
too  sensitive  to  edge  fine-structure  and  to  weak  gradients  to  be  useful  in  this  application.  The  edge  algorithm  pro¬ 
posed  by  Holyer  and  Peckinpaugh  [8]  is  based  on  the  cluster  shade  texture  measure,  which  is  derived  from  the 
gray  level  co-occurrence  (GLC)  matrix.  The  authors  [8]  have  suggested  that  the  edge  detection  technique  based 
on  the  GLC  matrix  can  be  effectively  used  in  automated  detection  of  mesoscale  features.  The  (/  ,j )  th  element  of 
the  GLC  matrix  ,  P(i  J  \  Ax  ,A y ),  is  the  relative  frequency  with  which  two  image  elements,  separated  by  distance 
(Ax  ,Ay ),  occur  in  the  image,  one  with  intensity  level  i  and  the  other  with  intensity  level  j.  The  elements  of  the 
GLC  matrix  could  be  combined  in  many  different  ways  to  give  a  single  numerical  value  that  would  be  a  measure 
of  the  edges  present  in  the  image.  Holyer  and  Peckinpaugh  [7]  have  used  a  cluster  shade  function  which  is  found 
to  be  very  effective  in  the  edge  detection  process.  The  new  edge  algorithm  computes  the  cluster  shade  function  at 
each  pixel.  Then  the  edges  are  detected  by  finding  the  significant  zero  crossings  in  the  cluster  shade  image.  The 
advantages  of  this  new  edge  algorithm  over  the  conventional  derivative-based  techniques  are  discussed  in  [7].  It  is 
known  that  using  large  windows  in  derivative-based  edge  detector  algorithms  results  in  poor  smoothing.  This 
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problem  is  circumvented  in  the  new  algorithm.  Because  edges  are  detected  by  finding  zero  crossings,  precisely 
positioned  lines  result,  even  it'  the  GLC  matrix  is  calculated  using  a  larger  window.  So,  the  desired  edge  detection 
characteristic  of  retaining  sharp  edges  while  eliminating  edge  detail  is  achieved  with  the  help  of  the  new  algo¬ 
rithm.  As  an  input  to  our  feature  labeling  algorithm,  we  used  the  image  output  generated  by  cluster  shade  algo- 
rithm.  with  a  window  size  of  16x16  pixels  and  zero  crossing  threshold  of  50.  The  edge  magnitudes  obtained  from 
tins  new  edge  detector  algorithm  are  used  as  an  input  to  the  feature  labeling  algorithm.  In  particular  the  edge  mag¬ 
nitudes  are  used  to  evaluate  the  a  priori  probability  values.  It  can  be  seen  from  figure  4.  that  only  the  edge  pixels 
participate  actively  in  the  labeling  process.  Also  the  speed  of  the  algorithm  that  computes  the  a  priori  probability 
values  is  greatly  improved  because  of  the  reduced  number  of  pixels  considered. 

In  the  second  stage,  the  iterative  updating  rule  is  implemented.  The  compatibility  coefficients  are  evaluated 
using  equation  (5.6).  The  coefficients  are  calculated  using  the  initial  class  probability  values  obtained  in  the  first 
stage.  These  are  fixed  during  the  update  process.  As  a  concluding  step  in  the  second  stage  of  the  implementation, 
the  iterative  updating  algorithm  is  implemented  using  the  equation  (5.9).  The  iterative  algorithm  terminates  when 
(p(~!  (x,:)  -pi(. r,r) )  <  £,  where  £  is  a  very  small  quantity.  Clearly  the  time  parameter  is  helpful  Li  reducing 
the  ambiguity,  especially  when  the  image  is  contaminated  wdth  noise. 


7.  PROTOTYPE  EXPERT  SYSTEM 

The  expert  system  described  here  is  based  on  an  embryonic  model  of  the  expected  movement  and  evolution 
of  warm-core  rings  (WCRs)  and  cold-core  rings  (CCRs).  The  informanon  in  the  knowledge  base  comes  from  a 
compilation  of  ruies  produced  by  a  survey  of  the  oceanographic  literature.  The  prototype  expert  system  is 
presendy  implemented  on  a  VAX  8300  in  a  combination  of  computer  languages.  The  rule  based  modules  are 
coded  in  OPS83.  Some  supporting  modules  that  provide  primarily  arithmetic  calculation  are  in  C.  The  module 
that  provides  the  graphics  support  is  in  FORTRAN. 


7.1  PROTOTYPE  EXPERT  SYSTEM  CAPABILITIES 

The  prototype  oceanographic  expert  system  depicts  expected  motions  and  size  changes  of  warm-core  rings 
(WCRs)  and  cold-core  rings  (CCRs),  and  their  interactions  with  each  other  and  with  the  Gulf  Stream,  in  different 
regions  of  the  northwest  Atlantic  Ocean.  The  domain  of  the  system  is  bounded  on  the  west  by  the  east  coast  of  the 
United  States  and  Canada.  There  are  no  explicit  limits  on  the  east,  north  and  south.  The  domain  is  divided  into 
nine  regions,  primarily  based  on  longitude.  The  rules  that  describe  the  expected  behavior  of  WCRs  and  CCRs  , 
and  Gulf  Stream  activity,  are  different  in  each  region. 

In  each  time  step,  each  ring  has  its  "status"  (i.e.,  position,  size)  updated.  That  operation  is  followed  by  a 
determination  of  the  ring’s  interaction  with  the  Gulf  Stream.  The  ring-Gulf  Stream  interaction  is  determined  by 
the  distance  and  direction  of  the  center  of  the  ring  to  the  nearest  point  on  the  Gulf  Stream  boundary.  Because  of 
the  ring-Gulf  Stream  interaction,  Gulf  Stream  evolution  is  also  simulated.  The  flow  chart  given  in  figure  5.  dep¬ 
icts  the  important  steps  involved  in  the  prototype  expert  system. 


8.  SUMMARY  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 


In  this  paper,  the  need  for  automatic  interpretation  of  oceanographic  images  is  emphasized.  The  advantage 
of  exploiting  the  contextual  information  in  feature  labeling  is  highlighted.  An  efficient  and  simple  technique  for 
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labeling  of  oceanic  features  is  described.  The  underlying  theoretical  framework  and  the  steps  involved  in  estimat¬ 
ing  the  probability  functions  are  explained  in  detail.  The  overall  architecture  of  the  proposed  image  interpretation 
system  is  presented.  Results  of  the  first  stage  of  the  labeling  technique  are  presented.  The  first  stage,  namely,  the 
estimation  of  a  priori  probabilities  for  the  reference  image,  is  complete  and  has  been  tested  on  real  oceanographic 
images.  An  expert  system  that  uses  the  knowledge  base  and  the  information  obtained  from  the  labeling  algorithm 
is  in  the  process  of  development. 

As  a  future  extension  to  the  present  work  we  propose  to  : 

(i)  investigate  the  possibility  of  implementing  a  parallel  relaxation  labeling  algorithm  to  speed  up  the  labeling 
process. 

(ii)  make  the  oceanographic  expert  system  "learn"  from  its  past  experience  in  analyzing  the  satellite  imagery 
data. 
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figure  1.  typical  IR  image  figure  2.  architecture  of  the  hybrid  system 
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figure  3.  hand  segmented  image 


figure  4.  results  of  the  labeling  algorithm  figure  5.  steps  of  prototype  expert  system 
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ABSTRACT 

Thermal  infrared  images  of  the  ocean  obtained  from 
satellite  sensors  are  widely  used  for  the  study  of  ocean  dynam¬ 
ics.  The  derivation  of  mesoscaie  ocean  information  from 
sateihte  data  depends  to  a  large  extent  on  the  correct  interpre¬ 
tation  of  infrared  oceanographic  images.  The  difficulty  of  the 
image  analysis  and  understanding  problem  for  oceanographic 
images  .$  uue  in  large  part  to  the  lack  of  precise  tnathemadcal 
descriptions  of  the  ocean  features,  coupled  with  the  time  vary¬ 
ing  nature  of  these  features  and  the  complication  that  the  view 
of  the  ocean  surface  is  typically  obscured  by  clouds,  some¬ 
times  almost  completely.  Towards  this  objective,  the  present 
paper  describes  a  technique  that  utilizes  a  non  linear  proba¬ 
bilistic  relaxation  method  for  the  oceanographic  feature  label¬ 
ing  problem.  A  unified  mathematical  framework  that  helps  in 
solving  the  problem  is  presented.  This  paper  highlights  the 
advantages  of  using  the  contextual  information  in  the  feature 
labeling  algorithm.  The  feature  labeling  technique  makes  use 
of  a  new,  efficient  edge  detection  algorithm  based  on  cluster 
shade  texture  measure.  This  new  algorithm  is  found  to  be 
more  suitable  for  labeling  the  mesoscaie  features  present  in 
the  oceanographic  satellite  images.  The  paper  presents  some 
important  results  of  the  series  of  experiments  conducted  at 
Remote  Sensing  Branch,  NORDA  on  the  NOAA  AVHRR 
imagery  data.  The  paper  concludes  with  a  motivation  for  using 
this  technique  to  build  an  oceanographic  expert  system. 

Key  words:  feature  labeling,  feature  extraction,  oceanic 
features,  edge  detection,  knowledge  based  systems,  relaxation, 
infrared  imagery. 


1.0  INTRODUCTION 

Satellite-bome  sensors  potentially  offer  many 
advantages  for  the  study  of  oceanic  processes.  They 
provide  global  synoptic  measurements  of  various 
oceanic  surface  properties,  in  contrast  to  the  local 
measurements,  possibly  at  a  range  of  depths,  pro¬ 
vided  by  conventional  oceanographic  measurement 
techniques.  Thermal  infrared  images  of  the  ocean 
obtained  from  satellite  sensors  are  widely  used  for 
the  study  of  ocean  dynamics.  Fig.l  shows  a  sample 
infrared  image  of  the  Gulf-Stream  obtained  from  the 
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Advanced  Very  High  Resolution  Radiometer 
(AVHRR)  aboard  the  NOAA-7  satellite.  Brightness 
in  this  infrared  image  is  inversely  proportional  to  the 
ocean  surface  temperature  (dark  areas  represent 
warmer  temperatures  and  light  areas  represent  colder 
temperatures).  Vortices  (areas  of  closed  circulation) 
within  this  turbulent  flow  pattern  are  called  eddies. 
The  Gulf  Stream  and  its  associated  eddies  are  exam¬ 
ples  of  mesoscaie  features  ("mesoscaie"  is  the  name 
commonly  applied  to  the  features  existing  on  spatial 
scales  of  the  order  of  50  to  300  km).  Mesoscaie 
features  are  important  to  the  study  of  ocean  dynam¬ 
ics,  the  fisheries  and  to  many  other  diverse  interests. 

Current  image  analysis  techniques  rely  on  human 
interpretation  of  the  satellite  imagery.  Human 
interpretation  is  obviously  varied  in  its  level  of 
expertise  and  is  highly  labor-intensive.  With  the  prol¬ 
iferation  of  high  volume  Advanced  Very  High  Reso¬ 
lution  Radiometer  image  applications  ,  it  becomes 
highly  desirable  for  certain  applications  to  move 
from  the  labor-intensive  manual  interpretation  of 
infrared  imagery  towards  a  capability  for  automated 
interpretation  of  these  images.  The  complete  automa¬ 
tion  of  the  oceanographic  image  interpretation  func¬ 
tion  is  probably  not  feasible,  but  one  can  begin  to 
address  certain  subsets  of  the  problem  with  the 
present-day  image  processing  and  artificial  intelli¬ 
gence  techniques.  TTiis  was  the  motivation  for  the 

work  reported  by  Lybanon  et  al„  [1]  (15],  in  the 
development  of  a  prototype  oceanographic  expert 
system. 

Several  previous  studies  have  addressed  the 
automation  of  the  analysis  of  the  infrared  satellite 
imagery  for  mesoscaie  features.  Gerson  and  Gabroski 
[2j  and  Gerson  et  al.,[3j  investigated  the  detection  of 
the  Gulf  Stream  in  infrared  images  from  the  Geosta¬ 
tionary  Operational  Environmental  Satellite  (GOES). 
Gerson  and  Gabroski  [2],  used  a  hierarchical 
approach  where  16x16  pixel  (128x128  km)  "frames" 
within  the  image  were  evaluated  for  the  possibility  of 
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containing  the  Gulf  Stream.  Frames  flagged  as  Gulf 
Stream  possibilities  were  then  further  evaluated  to 
determine  the  exact  location  of  the  Stream  within  the 
frame  by  looking  at  statistics  based  on  5x5  pixel 
"locai  neighborhood".  As  an  outgrowth  of  the  work 
reported  in  [21  and  [3],  Couiter  [4]  performed 
automated  feature  extraction  studies  using  the  higher 
resolution  Advanced  Very  High  Resolution  Radiom¬ 
eter  data.  Janowitz  [5]  studied  the  automatic  detec¬ 
tion  of  the  Gulf  Stream  eddies  using  the  Advanced 
Very  High  Resolution  Radiometer  data.  Nichol  (6) 
uses  a  region  adjacency  graph  to  define  spatial  rela¬ 
tionships  between  elementary  connected  regions  of 
constant  gray  level  called  atoms.  Eddy-like  structure 
is  then  identified  by  searching  the  graph  for  isolated 
atoms  of  high  temperature  that  are  enclosed  by  atoms 
of  lower  temperature  (for  the  case  of  warm  eddies). 
Although  satisfactory  emulation  of  human  extraction 
of  eddy  structure  is  claimed  for  this  method. 
NichoI[61  does  point  out  that  not  all  enclosed  uniform 
areas  identified  by  the  method  will  correspond  to  real 
ocean  structure. 


2.0  MOTIVATION  OF  THE  PRESENT  WORK 

In  this  section  we  present  the  motivation  behind 
the  present  work.  Our  primary  objective  is  to  build  a 
powerful  automatic  image  interpretation  system  for 
oceanographic  satellite  images.  In  order  to  make  this 
difficult  problem  tractable,  we  divide  the  problem 
into  two  parts  :  feature  labeling  problem  and  the 
development  of  an  expen  system.  In  this  paper  we 
focus  on  the  feature  labeling  problem.  It  is  dear  that 
the  performance  of  the  labeling  algorithm  depends 
heavily  on  the  low  level  image  processing  algo¬ 
rithms.  Particularly,  the  output  of  an  edge  detector 
algorithm  plays  a  major  role  in  the  feature  labeling 
process.  In  view  of  this,  a  new  efficient  edge  detector 
algorithm  proposed  by  Holyer  and  Peckinpaugh  [7]  is 
employed  in  our  feature  labeling  technique,  ft  is 
know-n  that  the  conventional  edge  derivative  opera¬ 
tors  are  very  sensitive  to  noise  and  are  not  suitable 
for  analyzing  oceanographic  satellite  images.  The 
new  edge  detector  algorithm  proposed  by  Holyer  and 
Peckinpaugh  [7]  is  based  on  the  gray  level  co¬ 
occurrence  matrix,  which  is  commonly  used  in  image 
texture  analysis.  This  algorithm  (71  is  found  to  exhi¬ 
bit  the  characteristics  of  fine  structure  rejection  while 
retaining  edge  sharpness.  In  this  paper  we  focus  on 
the  preliminary  results  of  a  feature  labeling  algorithm 
and  the  effect  of  the  duster  shade  edge  detector  algo¬ 
rithm  on  the  performance  of  the  labeling.  In  the  fol¬ 
lowing  paragraph  we  briefly  outline  the  various  steps 
taken  to  solve  the  second  sub  problem. 


The  objective  of  an  expert  system  for  oceano¬ 
graphic  system  is  to  correctly  interpret  the  dynamics 
of  the  ocean  process  with  minimal  human  interaction. 
Towards  this  objective,  the  development  of  a  power¬ 
ful  oceanographic  expert  system  is  in  the  process  of 
development  at  Remote  Sensing  Branch,  N'ORDA. 
We  have  already  developed  a  prototype  expert  sys¬ 
tem  for  Gulf  Stream  regional  dynamics  [1],  Tne 
current  activities  focus  on  building  an  expert  system 
that  makes  use  of  the  information  about  the  position 
cf  mesoscale  features  obtained  from  the  relaxation 
labeling  scheme.  One  element  of  an  expert  system  is 
a  database  of  knowledge  about  the  subject  matter, 
with  the  knowledge  represented  in  a  form  suitable  fc. 
manipulation  by  the  "inference  engine"  of  me  expert 
system  (i.e.,  the  logical  and  heuristic  procedures  for 
solving  problems  in  the  pmblem  domain).  With  the 
help  of  the  knowledge  gained  from  discussions  at 
NORDA  and  from  the  literature  in  oceanography,  a 
knowledge  base  about  the  nature  of  the  mesoscale 
features  (occurrence,  mean  life  time,  movement,  etc.) 
was  built  [1].  Presently,  we  are  in  the  process  of 
designing  an  expert  system  which  makes  use  of  the 
positional  information  about  the  featurestobtained 
from  the  relaxation  labeling)  and  the  knowledge  base 
of  mesoscale  information. 

This  paper  presents  a  relaxation  labeling  scheme 
to  label  features  in  thermal  infrared  images  obtained 
from  satellite.  The  technique  presented  in  this  paper 
exploits  the  advantages  of  using  the  contextual  infor¬ 
mation  in  the  labeling  algorithm.  Tne  remainder  of 
the  paper  is  organized  as  follows:  Section  3.0  intro¬ 
duces  the  probabilistic  relaxation  scheme  and  briefly 
discusses  the  application  of  this  scheme  to  the 
oceanographic  labeling  problem.  In  section  4.0  we 
develop  a  mathematical  frame  work  which  is  neces¬ 
sary  for  our  labeling  problem.  Section  5  0  discusses 
the  performance  of  our  technique  and  describes  the 
steps  that  are  involved  in  implementing  this  tech¬ 
nique.  Section  6.0  gives  a  brief  description  on  the 
implementation  of  this  technique  on  the  image 
analysis  system  at  NORDA.  Section  7.0  concludes 
with  ail  the  important  features  of  this  technique  and 
future  extensions  to  the  available  scheme. 


3.0  RELAXATION  PROCESS 

An  important  research  area  in  image  analysis  and 
image  interpretation  technology  is  the  development 
of  methods  that  blend  contextual  information  with 
the  conventional  image  processing  algorithms.  The 
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literature  survey  clearly  indicates  that  such  a  hybrid 
approach  yields  good  results.  Relaxation  labeling  is 
one  such  process  that  approaches  this  problem. 
Relaxation  labeling  has  been  applied  to  a  variety  of 
image  processing  problems  eg.,  linear  feature 
enhancement  [8],  edge  enhancement  [9j,  image 
enhancement  [10],  pixel  classification  [11,12]. 
Recently  a  survey  article  by  Kitten  and  Illingworth 
[13]  on  relaxation  labeling  highlights  the  importance 
of  this  area  of  research.  The  survey  [131  also  points 
cut  the  advantages  and  possible  applications  of  relax- 
ation  methods.  More  importantly,  the  relxxation 
labeling  approach  was  elegantly  described  by  Rosen- 
feld  et  al.  [14]  who  investigated  the  problem  of  label¬ 
ing  the  sides  of  a  triangle  and  proposed  a  set  of 
schemes  to  solve  the  problem.  The  paper  [14]  con¬ 
cluded  with  the  result  that  the  non  linear  probabilistic 
relxxation  schemes  yield  better  results  than  the  eth¬ 
ers.  The  labeling  algorithm  presented  in  this  paper  is 
based  on  the  non  linear  probabilistic  relxxation  tech¬ 
nique. 


The  goal  of  the  relaxation  process  is  to  reduce 
the  uncertainty  (and  improve  the  consistency)  in  the 
assignment  of  one  of  the  labels  to  each  object  in  a  set 
of  related  objects.  In  the  Oceanographic  feature 
classification  problem,  the  classes  are  the  various 
Oceanographic  features  namely  North  and  South  wall 
of  Gulf  Stream,  cold  eddy,  warm  eddy,  shelf  front 
and  coastal  boundary.  Refer  to  fig.2  to  identify  the 
position  of  these  oceanic  features  in  a  typical  image. 
The  objects  are  the  individual  pixels  in  a  set  of 
registered  multi-temporal  images.  The  uncertainty 
could  be  due  to  the  cloud  cover  or  the  overlap  of  the 
features  that  are  mentioned  above,  features  not 
belonging  to  one  of  the  classes,  noise  in  the  image,  or 
other  factors.  In  this  paper,  we  are  attempting  to  label 
mesoscale  features,  but  the  ocean  exhibits  variability 
on  ail  spatial  scales.  Thermal  structure  on  scales 
smaller  than  mesoscale  will  interfere  with  the  mesos¬ 
cale  feature  labeling  process.  The  underlying 

mathematical  ffiamewoik  necessary  for  the  relaxation 
labeling  method  is  described  in  the  next  section  of 
the  paper. 


4.0  MATHEMATICAL  FRAMEWORK 

Let  A  =  }  be  the  set  of  possible  labels  that 

may  be  assigned  to  each  pixel  x  in  the  IR  image. 
Also  we  let  p{{x,t)  denote  the  probability  that  the 
pixel  at  x(i  j)  and  time  instant  t,  belongs  to  the  object 
X  after  k  iterations  of  the  relaxation  algorithm.  Note 


that  the  probabilities  are  functions  of  time  unlike  the 
conventional  pixel  relaxation  labeling  schemes  where 
die  probability  is  a  function  of  position  alone.  This 
allows  the  relaxation  labeiing  algorithm  to  utilize 
temporal  continuity  to  reduce  the  ambiguity  in  label¬ 
ing.  The  ambiguity  may  arise  due  to  noise  (cloud 
cover,  for  example). 


Tnere  are  two  steps  in  executing  the  probabilistic 
relaxation  algorithm.  At  the  first  step,  a  priori  proba¬ 
bilities  are  evaluated  with  the  help  of  ground  truth 
data  and  /  or  a  previous  but  recent  mesoscale 
analysis.  In  the  second  step,  these  a  priori  probabili¬ 
ties  are  iteratively  updated  (relaxation)  until  a  con¬ 
sistent  labeling  is  reached.  We  now  discuss  these  two 
steps  in  detail. 


Step  1:  Estimating  the  a  priori  probabilities 


Let  pxti.i)  denote  the  a  priori  value,  that  is,  the  pro¬ 
bability  that  pixel  x(i  j)  at  time  t  belongs  to  the  object 
X,  at  the  zeroth  iteration.  The  Bayesian  probability 
equation  is  used  to  evaluate  this  value.  The  equation 
(4.1)  is  used  to  calculate p£(x.:). 


P 


x  u  )  = 


p(x  j  IX)  PfX) _ 

Z\  P(x.t  IX)  P(X) 


l-  !) 


where  p(x.t  IX)  denote  the  conditional  density  function 
and  P(X)  the  probability  of  occurrence  of  the  object  X. 

To  evaluate  the  conditional  density  function  puj  iX), 
a  set  of  parameters  is  measured  at  the  pixel  x(i  j).  Let 
X  denote  the  parameter  vector.  The  following  param¬ 
eters  are  used  to  form  the  vector  X: 

(1)  distance  of  the  pixel  x(ij)  from  the  origin,  both 
the  magnitude  and  direction. 

(2)  gray  scale  intensity  value  at  the  pixel  x(i  j). 

(3)  the  edge  magnitude  (Section  5.1  presents  the 
chosen  edge  operator  algorithm). 

For  each  object,  the  mean  vector  m.  and  the  covari¬ 
ance  matrix  lx  are  computed.  Also  it  is  assumed  that 
the  conditional  density  function  follows  a  normal  dis¬ 
tribution.  Hence  the  conditional  density  function 
p(x,t  iX)  is  evaluated  using  equation  (4.2). 


P(XJ  U)  =  (2jiI£xI)-,W  -y  (X-Ux)'Ix'’(X-4x)r  (4-2> 
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To  compute  PCo  .  relative  areas  of  the  objects  are  cov(XX) 

considered.  The  number  of  pixels  La  the  object  X  be  ru(x.y)  =  —  - 

,  ,  4  ^ .  u(  A;  C(  A. ; 

v,.  Then  .»  X>  can  be  calculated  using  equation  (-».o). 

where  cov(k.X)  denotes  the  covariance  of  two  events, 
namely,  the  pixel  x  is  assigned  the  label  X  and  the 
p  cM  =  — (4.3)  pLxel  y  is  assigned  the  label  X  and  is  given  as: 


Step  2:  Iterative  updating  algorithm 

We  now  discuss  the  probability  updating  rule.  The 
new-  estimate  of  the  probability  of  X  at  x(ij)  is  given 

j . 


? 


I  1  ( 1  -ry  1  t'x  )) 


(I  -aO?xU/X4.4) 


cov  {XX )  ~  p  (X<&  X )  -  p  (X)p  (X )  (4.7) 


The  joint  probability  p(X&X)  can  be  calculated  using 
the  position  vectors  of  the  pixels  x  and  v.  Let  m  and 
m  be  the  mean  position  vectors  of  the  objects  X  and 
X’  respectively.  Then  the  joint  probability  can  be 
computed  using  simple  functions  as  given  by  equa¬ 
tions  (4.8a)  and  (4.Sb). 


where  el  :  is  called  the  update  factor  and  ax  is  called 
tempera!  weighting  function. 

Hie  method  of  estimating  these  factors  are  illustrated 
below.  The  updating  factor  for  the  estimate  plu./)  at 
the  >th  iteration  is  given  by  equation  (4.5). 


p(X  &  X)  =  — ~  if  X  *  X  (4.8a) 

m  -  m 

pCX  &  X)  -  -  — — if  X  =  X  (4.8b) 

.Z  ?  iT» 


q  {■ x  )  =  -  Z,  Sv  r  u<x  .y  )p  i  O'  >  (4.5) 


where  m  is  the  number  of  objects.  In  this  equation, 
'■«<•*■/>  denote  compatibility  coefficients.  These 
coefficients  are  computed  as  in  (13,16).  According  to 
the  relaxation  scheme,  r^ix.y)  is  a  measure  of  the 
probabilistic  compatibility  between  label  X  on  point  x 
and  label  X  on  point  y,  and  has  the  following  charac¬ 
teristics: 

(  I )  If  X  on  x  frequently  co-occurs  with  X'  on  y,  then 
rxx<i.y)  >  0  ,  and  if  they  always  co-occur,  then 
r\x<x-y)  =  1- 

(2)  If  X  on  x  rarely  co-occurs  with  X  on  y,  then 
rurix.y)  <  0  -  4™!  if  they  never  co-occur,  then 

rix <*-y)  ~  -1- 

(3)  If  X  on  x  occurs  independently  of  X  on  y,  then 
rua,y)  =  0. 

These  compatibility  coefficients  are  computed  using 
the  equation  (4.6). 


The  new  estimate  of  the  probability  of  X  at  x(i\j)  is  : 


pl*\x.t)  =<X\ 


p{{.Xj)<\*ql<x)) 

I  mi(x,;)  (Hihl) 


+  (t-a\)p\(.x  /  f4.9) 


where  pxtxS)  denotes  the  probability  value  at  the 
pixel  x(ivj)  at  time  instant t’  <  t  and  ax  denotes  the 
temporal  weighting  function  for  the  object  X.  The 
function  ax  determines  the  weight  that  is  given  to  the 
current  probability  values  and  the  factor  (l-cu)  associ¬ 
ates  a  weight  to  the  probability  values  calculated  at 
previous  time  instant  r'. 

To  estimate  the  new  probability  value,  either  (a) 
the  probability  value  at  the  previous  time  frame  t’  or 
(b)  a  set  of  probability  values  at  time  frames  t’ ,  t'  , 

etc.,  can  be  used.  In  the  present  analysis  and  the 
implementation,  the  method  (a)  is  chosen.  Also  it  is 
proposed  to  implement  the  method  (b)  and  compare 
improvements  (if  any)  in  the  feature  labeling.  The 
temporal  weighting  function  can  be  determined  with 


the  he of  the  previous  analysis  of  the  images.  It  is 
understood  that  this  temporal  weighting  function  can 
not  be  constant  in  the  real  time  situation.  For  example 
in  some  cases,  the  performance  of  the  labeling  may 
be  less  satisfactory  and  correspondingly  the  temporal 
sighting  function  should  be  given  a  relatively 
smaller  value.  Hence  the  selection  of  this  function 
itself  may  be  a  problem  in  some  cases.  One  may  use 
different  values  for  this  function  and  the  performance 
of  the  labeling  can  be  compared.  The  iterative  updat¬ 
ing  is  terminated  when  the  difference  between  the 
probability  values  at  the  k  th  and  the  (k  +  I)  th 
iterations  is  very  small  (say  less  than  0.1  Tc). 

5.0  DISCUSSION  ON  THE  TECHNIQUE 

Tne  implementation  of  the  above  mentioned 
techr.iuue  is  carried  out  in  two  stases.  The  flow  chart 
shown  in  fig. 4  depicts  the  various  steps  involved  in 
implementing  the  technique. 


At  the  first  stage,  the  a  priori  probabilities  are 
estimated  using  a  manually  prepared  mrsoscale 
analysis  from  a  time  period  of  five  days  prior  to  the 
test  image. 

The  objects  present  in  the  oceanographic  IR 
image  are  identified  with  the  help  of  the  previous 
analysis.  The  shapes  (or  boundaries)  of  these  objects 
are  then  determined  and  represented  by  a  regular 
polygon.  The  parameter  values  are  then  computed  for 
all  the  pixels  inside  these  polygons.  A  reasonably 
accurate  algorithm  which  finds  out  whether  a  pixel  is 
inside  a  polygon  or  not  is  implemented. 

From  the  parameter  vector  X  for  a  pixel  x(ij), 
the  mean  vector  p*  is  calculated.  Also  for  each  object 
X  ,  the  covariance  matrix  I*  is  computed.  The  equa¬ 
tion  (4.2)  is  used  to  compute  the  conditional  density 
function  p(x,tW).  Finally  the  initial  probability 
p°(x,t )  is  computed  using  the  equation  (4.1).  The 
output  of  the  first  stage  is  shown  in  fig.3.  The  output 
image  is  split  into  six  parts  to  illustrate  the  labeling. 
Fig.3. 1  shows  the  output  of  the  labeling  module  with 
all  the  five  features  labeled.  Refer  to  fig. 2  to  compare 
the  labeling  with  the  hand  segmented  image.  To  illus¬ 
trate  the  labeling  process  five  different  features  are 
considered.  They  are  :  two  parts  of  the  north  wall, 
south  wall  and  two  warm  eddies.  The  labeling  of 
these  five  features  is  shown  in  figures  3.2  to  3.6.  The 


individual  pixels  in  the  oceanographic  image  arc 
labeled  depending  on  the  initial  probability  values 
obtained  in  the  first  stage. 

5.1  DISCUSSION  ON  CLUSTER  SHADE  EDGE 

ALGORITHM 

In  this  section,  we  discuss  the  features  of  the 
edge  detection  algorithm  proposed  by  Holver  and 
Peckinpaugh  [7 J.  Tne  motivation  behind  the  develop¬ 
ment  of  such  a  new  edge  detector  algorithm  is  to  aid 
the  analysis  of  oceanographic  satellite  images.  The 
popular  derivative-based  edge  operators  viz  Sobel’s 
operator  [19]  are  shown  to  be  too  sensitive  to  edge 
fine-structure  and  to  weak  gradients  to  be  useful  in 
this  application.  Tne  edge  algorithm  proposed  by 
Holver  and  Peckinpaugh  [7]  is  based  on  the  cluster 
shade  texture  measure,  which  is  derived  from  the 
gray  level  co-occurrence  matrix  (GLC).  Tne  authors 
[7]  have  suggested  that  the  edge  detection  technique 
based  on  the  GLC  matrix  can  be  effective!',  used  in 
automated  detection  of  mesoscale  features.  The  U.j ) 
th  element  of  the  GLC  matrix  ,  P(i.j  !4u.4y  s  is  the 
relative  frequency  wuth  which  two  image  elements, 
separated  by  distance  (At.Av)  occur  in  the  image,  one 
with  intensity  level  i  and  the  other  with  intensity  level 
j.  The  elements  of  the  GLC  matrix  could  be  com¬ 
bined  in  many  different  ways  to  give  a  single  numeri¬ 
cal  value  that  would  be  a  measure  of  the  edges 
present  in  the  image.  Holver  and  Peckinpaugh  [7] 
have  used  a  cluster  shade  function  which  is  found  to 
be  very  effective  in  the  edge  detection  process.  The 
new  edge  algorithm  computes  the  cluster  shade  func¬ 
tion  at  each  pixel.  Then  the  edges  are  detected  by 
finding  the  significant  zero  crossings  in  the  cluster 
shade  image.  The  advantages  of  this  new  edge  algo¬ 
rithm  over  the  conventional  derivative-based  tech¬ 
niques  are  discussed  in  [7].  It  is  known  that  using 
large  windows  in  derivative-based  edge  detector 
algorithms  results  in  poor  smoothing.  This  problem  is 
circumvented  in  the  new  algorithm.  Because  edges 
are  detected  by  finding  zero  crossings,  precisely  posi¬ 
tioned  lines  result,  even  if  the  GLC  matrix  is  calcu¬ 
lated  using  a  larger  window.  So,  the  desired  edge 
detection  characteristics  of  retaining  sharp  edges 
while  eliminating  edge  detail  is  achieved  with  the 
help  of  the  new  algorithm.  As  an  input  to  our  feature 
labeling  algorithm,  we  used  the  image  output 
generated  by  cluster  shade  algorithm,  with  a  window 
size  of  16x16  pixels  and  zero  crossing  threshold  of 
50.  The  edge  magnitudes  obtained  from  this  new 
edge  detector  algorithm  is  used  as  an  input  to  the 
feature  labeling  algorithm.  In  particular  the  edge 
magnitudes  are  used  to  evaluate  the  a  priori  probabx  - 
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icy  values  Ic  can  be  seen  from  fig  3.  that  the  only  the 
eJge  pixels  participate  actively  in  che  labeling  pro¬ 
cess  .  Also  the  speed  of  the  algorithm  that  computes 
the  a  priori  probability  values  is  greatly  improved 
because  of  the  reduced  number  of  pixels  for  con¬ 
sideration. 

In  the  second  stage,  the  iterative  updating  rule  is 
implemented.  The  compatibility  coefficients  are 
evaluated  using  the  equation  (4.6).  The  coefficients 
are  calculated  using  the  initial  class  probability 
values  obtained  in  che  first  stage.  These  are  fixed  dur¬ 
ing  the  update  process.  As  a  concluding  step  in  the 
second  stage  of  the  implementation,  the  iterative 
updating  algorithm  is  implemented  using  the  equa¬ 
tion  (4.9).  The  iterative  algorithm  terminates  when 
(x .:)  -p{  ix.:) )  <  e  where  e  is  a  very  small 
quantity.  Clearly  the  time  parameter  is  helpful  in 
reducing  the  ambiguity,  especially  when  the  image  is 
contaminated  with  noise. 


6.0  IMPLEMENTATION 

In  this  section,  we  present  the  experimental 
results  of  the  first  stage  of  the  relaxation  labeling 
technique  :  estimating  the  a  priori  probability  values 
of  a  reference  image.  The  software  that  computes  the 
a  priori  probabilities  is  developed  on  VAX  8300  sys¬ 
tem  running  the  VMS  operating  system  installed  at 
the  Remote  Sensing  Branch,  NORDA.  We  used  the 
Interactive  Digital  Satellite  Image  Processing  System 
consisting  of  International  Imaging  System’s  S600 
software  with  N'ORDA  extensions  to  process  the 
oceanographic  satellite  images.  The  software 
modules  are  developed  in  ’C’  language.  The  posi¬ 
tional  information  about  the  oceanographic  features 
present  in  the  image  can  be  provided  either  manually 
or  from  the  previous  analysis  of  the  satellite  imagery 
data.  The  module  first  computes  the  mean  vectors 
and  covariance  matrices  for  all  the  objects  present  in 
the  reference  image.  The  second  part  of  the  software 
module  computes  the  a  priori  probability  values  as 
given  by  equation  (4.1). 

To  illustrate  the  performance  of  the  labeling 
algorithm,  an  oceanographic  infrared  image  with 
some  of  the  typical  oceanic  features  is  taken  (refer  to 
fig.l).  The  objects  that  are  selected  for  labeling  are 
north  wall,  south  wall  and  warm  eddies.  The  a  priori 
knowledge  about  the  objects,  namely  the  position  is 
taken  from  the  prior  analysis  shown  in  fig. 3  and 
given  as  an  input  to  the  labeling  module.  The  a  priori 
probability  values  are  then  computed  for  all  the  pix¬ 


els  present  in  the  image.  To  highlight  the  perfor¬ 
mance  of  this  initial  labeling  algorithm,  a  decision 
making  algorithm  based  on  these  a  priori  probability 
values  is  implemented.  A  pixel  that  has  a  probability 
value  above  0.8  for  an  object  is  given  the  label  for 
that  object.  The  decision  algorithm  selects  the  pixels 
in  the  neighborhood  of  the  objects  to  illustrate  the 
results  obtained.  The  edges  of  the  labeled  objects  are 
dilated  for  clear  identification  of  the  labels.  Tne 
covariance  matrices  for  the  objects  are  shown  in 
fig.5. 


7.0  SUMMARY  AND  RECOMMENDATIONS 
FOR  FUTURE  WORK 

In  this  paper,  the  need  for  automatic  interpreta¬ 
tion  of  oceanographic  images  is  emphasized.  Tne 
advantage  of  exploiting  the  contextual  information  in 
feature  labeling  is  highlighted.  An  efficient  and  sim¬ 
ple  technique  for  labeling  of  oceanic  features  is 
described.  The  underlying  theoretical  framework  and 
functions  are  explained  in  detail.  Results  of  the  first 
srage  of  the  labeling  technique  are  presented.  The 
results  and  the  performance  analysis  of  the  proceed¬ 
ing  stages  of  the  technique  will  appear  in  forthcom¬ 
ing  issues.  The  first  stage,  namely,  the  estimation  of  a 
priori  probabilities  for  the  reference  image,  is  com¬ 
plete  and  has  been  tested  on  real  oceanographic 
images.  An  expert  system  that  uses  the  knowledge 
base  and  the  information  obtained  from  the  labeling 
algorithm  is  in  the  process  of  development. 

As  a  future  extension  to  the  present  work  we 
propose  to  : 

(i)  investigate  the  possibility  of  implementing  a 
parallel  relaxation  labeling  algorithm  to  speed  up 
the  labeling  process. 

(ii)  make  the  oceanographic  expert  system  "learn" 
from  its  past  experience  in  analyzing  the  satellite 
imagery  data. 

(iii)  incorporate  additional  information  about  the 
mesoscale  features  (size,  physical  properties, 
etc.)  during  the  labeling  process  to  reduce  the 
ambiguity  in  labeling  the  features. 
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fig-3.2  -  3.6  individual  features  labeled 
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